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Abstract. The analysis of high-energy air shower data allows one to study the 
proton-air cross section at energies beyond the reach of fixed target and collider 
experiments. The mean depth of the first interaction point and its fluctuations are a 
measure of the proton-air particle production cross section. Since the first interaction 
point in air cannot be measured directly, various methods have been developed in 
the past to estimate the depth of the first interaction from air shower observables 
in combination with simulations. As the simulations depend on assumptions made 
for hadronic particle production at energies and phase space regions not accessible in 
accelerator experiments, the derived cross sections are subject to significant systematic 
uncertainties. The focus of this work is the development of an improved analysis 
technique that allows a significant reduction of the model dependence of the derived 
cross section at very high energy. Performing a detailed Monte Carlo study of the 
potential and the limitations of different measurement methods, we quantify the 
dependence of the measured cross section on the used hadronic interaction model. 
Based on these results, a general improvement to the analysis methods is proposed by 
introducing the actually derived cross section already in the simulation of reference 
showers. The reduction of the model dependence is demonstrated for one of the 
measurement methods. 
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1. Introduction 

The natural beam of cosmic ray particles extends to energies far beyond the reach of any 
Earth-based particle accelerator. Therefore cosmic ray data provide a unique window 
to study hadronic interaction phenomena at energies up to several Joules per particle, 
corresponding to an equivalent center-of-mass energy of up to 450 TeV. 

On the other hand, direct observation of the first interaction of ultra-high energy 
cosmic rays in the upper atmosphere is impossible due to the very low flux of these 
particles. Only the cascades of secondary particles, called extensive air showers (EAS), 
can be measured with arrays of particle detectors or optical telescopes. To obtain 
information on the first interactions in an air shower it is necessary to link the measured 
air shower characteristics to that of high energy particle production in the shower. This 
can be done with detailed Monte Carlo simulations of the shower evolution and the 
corresponding shower observables, but inevitably causes a dependence of the results on 
hadronic interaction models needed for the shower description. 

The total particle production cross section is one of the most fundamental quantities 
that characterizes hadronic interactions. Considering proton-induced air showers of the 
same primary energy, the depth of the first interaction point, Xi, is distributed according 
to 



dXi Ap_air 

where Ap_air is the interaction mean free path of protons in air. The mean depth of the 
first interaction point and its shower-to-shower fluctuations are directly linked to the 
size of the proton-air cross section by 

L'p-air — , ) {■^J 

^p— air 

with the mean target mass of air being (mair) ~ 14.45 rrip = 24160 mbg/cm^ [1]. It 
is, therefore, not surprising that there is a long history of attempts to infer this cross 
section from high-energy cosmic ray data [2-15]. 

A compilation of published proton-air cross section measurements and predictions of 
ultra-high energy hadronic interaction models is shown in Fig. [H All data above 100 TeV 
are based on the analysis of EAS data [7-13]. Also the ARGO-YBJ measurements 
around 10 TeV are originating from a high altitude EAS array [14]. The data at lower 
energies stem from unaccompanied hadron analyses [2-6] . 

All cosmic ray measurements of the proton-air cross section are only sensitive to 
the particle production cross section [9,22]. In addition, interactions with insignificant 
particle production have no measurable impact on cosmic ray observations. This is the 
case for air shower based techniques as well as for the unaccompanied hadron method. 

To draw conclusions on the energy of the primary particle and its first ultra-high 
energy interactions, all relevant processes involved from the primary cosmic ray particle 
entering the atmosphere up to the measurement of the EAS observable need to be 
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Figure 1. Compilation of proton-air production cross sections from cosmic ray 
measurements [2-14]. The data are compared to model predictions [16-21]. 



modeled as precisely as possible. Sophisticated EAS Monte Carlo simulation programs 
are available for this task. This results in a highly indirect analysis procedure. 

In this article we will perform a detailed Monte Carlo study of different methods 
to derive the proton-air cross section from air shower data at an energy of ~ 10^^ eV. 
We will assume protons as cosmic ray particles in the energy range considered here. 
Although there exist theoretical models [23-27] and also experimental indications [28-30] 
for this flux being indeed dominated by protons, other elements in the primary flux have 
also to be taken into account. This will be done in a forthcoming article that will address 
specifically this topic [31]. 

Based on the results of the Monte Carlo study of existing analysis methods a 
general improvement is proposed and explicitly applied to one of the methods. By 
accounting for the actually measured cross section already in the simulation of showers, 
the systematic uncertainty due to the limited knowledge of hadronic interactions is 
significantly reduced. The performance of the improved method is thoroughly tested for 
the application to high quality data of the depth of the shower maximum. Sources of 
systematic uncertainties of the resulting cross section are discussed. It is shown that the 
dependence on the hadronic interaction model, the most important source of systematic 
uncertainties, can be significantly reduced by incorporating the measured cross section 
in a consistent way in the shower simulation. 

The analyses of this work are done at an energy of 10^^ eV, at which high statistics 
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Figure 2. Definition of variables to characterize EAS longitudinal shower 
development. Zero slant depth is where the cosmic ray particle enters the atmosphere. 
The first interaction occurs at Xi. At Xmax the shower reaches its maximum particle 
number A'niax- After the maximum the shower is attenuated over AX2 before it reaches 
the particle number A'^gtage at the slant depth of Xobs- 

data are available or expected from HiRes [32], the Pierre Auger Observatory [33], and 
Telescope Array [34], but the results are qualitatively also valid at other energies. The 
expected reduction of the model dependence will be smaller at energies where data from 
accelerators are available, i.e. 10^^ eV and below. 

2. Relation between air shower observables and depth of the first 
interaction point 

Interactions over many decades in energy are occurring during EAS development. In the 
startup phase of the shower, relatively few ultra-high energy hadronic interactions are 
distributing the energy of the primary cosmic ray particle to a quickly growing number 
of secondary particles. The stochastic nature of these initial interactions is the source 
of strong fluctuations of EAS initiated by identical cosmic ray particles (shower-to- 
shower fluctuations). A significant part of secondary particles decay to electromagnetic 
particles and lead to the development of an electromagnetic cascade that, already after 
a few interactions, contains most of the total EAS energy and constitutes by far the 
largest fraction of the particles. The interactions of particles in the e.m. cascade are 
theoretically well understood. Also the large number of these interactions levels out 
additional large-scale fluctuation^. Thus, the development of an air shower can be 
characterized by two main stages: 

(i) Startup phase, consisting of few initial hadronic interactions at ultra-high 
energies. These interaction are the main source of fluctuations. 

II An exception are electromagnetic showers of i? > 10^^ eV that, by chance, happen to start very deep 
in the atmosphere, for which the Landau-Pomeranchuk-Migdal effects leads to very large shower-to- 
shower fluctuations [35-37]. 
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(ii) Cascade phase, in which very numerous interactions of the bulk of the shower 
particles at intermediate energies take place. No significant large-scale fluctuations 
are expected from this part. 

A clear boundary between the two stages cannot be drawn. The transition is seamless 
and by itself subject to strong shower-to-shower fluctuations. Only treatment with full 
air shower Monte Carlo simulation programs can fully account for all fluctuations. 

To discuss the relation of the depth of the first interaction point to air shower 
observables it is useful to introduce a simple longitudinal model of air shower 
development that connects the proton-air cross section to the observables in a 
transparent way. The naming conventions for the model are given in Fig. O We 
distinguish between two types of air shower observation: by ground based detector 
arrays and the direct observation of longitudinal shower profiles by telescope detectors. 
The relevant air shower observables used in proton-air cross section analyses are the 
position of the shower maximum, X^ax? or the combination of the number of electrons 
A'^e and muons iV^ at the atmospheric slant depth of the detector Xobs- 

In the following, the correlation of these EAS observables to the characteristics of 
the ultra-high energy interactions is studied with the one dimensional air shower Monte 
Carlo program CONEX v2r2 [38]. 

The particle numbers are defined as the total number of electrons above 1 MeV and 
muons above 1 GeV. This corresponds to typical quantities observed in air shower arrays, 
however the exact definition of observables depends very much on the experimental 
setup and varies strongly from experiment to experiment. By using the total number of 
particles for our study we are focusing on general air shower properties. 

The shower maximum Xmax is the slant depth of the maximum energy deposit of 
the shower in the atmosphere. This definition matches the shower maximum derived 
from the fluorescence light profile of showers and coincides with that of the particle 
number within ~ 3gcm~^. 

2.1. Arrays of particle detectors 

Using ground based air shower arrays, one can estimate the proton-air cross section by 
measuring the frequency of air showers of the same energy and stage of their development 
at different atmospheric depths. By selecting events of the same energy but different 
directions of incidence, the point of the first interaction has to vary with the angle, 
in order to observe the shower at the same stage of development. The selection of 
showers of constant energy and stage can be done only approximately and depends 
on the particular detector setup, but the typical requirement is a constant (Ae, A"^) at 
observation level. Examples of measurements of this type are [7-10,14]. By requiring 
a given number of muons at the detector level does, in first approximation, select EAS 
of the same primary energy, because the attenuation of muons is small. However, also 
muons are slowly attenuated in the atmosphere. To correct for that a constant intensity 
selection can be applied [39]. Showers with identical primary energy at the same stage 
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Figure 3. Distribution of the depth distance from the point of the first interaction up 
to the shower maximum for proton primaries at 10^^ eV. The predictions of different 
models are compared. 



of their shower development are assumed to yield the same number of electrons since, 
after the shower maximum, the electromagnetic shower attenuation is approximately 
universal [40-45]. 

The number of showers A'^sh selected per area dA and time dt, requiring a constant 
{Ne, N^) at the atmospheric depth Xobs can be written as 

— — — = dndE dXi dAXi dAX2 — 

dA^edAT^dXobsdAdt J Aint 

X PaxAE) Pax,{E,X,,AXi) Pn^{E,XuAXu. 

X 5(Xobs - Xi - AXi - AX2) 

^ dNcR 

^ dndEdAdt ' 

where dNcB./{dQdE dAdt) is the flux of cosmic ray particles. The distribution 
Paxi = dpi/dAXi describes shower-to-shower fluctuation of AXi, see Fig. |3l The two- 
dimensional probability density Pax2 = dp2/(dAX2 dXe) is the frequency of showers of 
given energy E, Xi, and AXi to have a shower size at a depth distance AX2 from the 
shower maximum. It is required that the electron number after the shower maximum 
is attenuated to N^, while Xobs = Xi + AXi + AX2. The selection of showers by 
energy according to their muon number is reflected by the energy dependent probability 
distribution P^^ = dp^/dN^. 

It is important to notice that the distribution Paxi does not depend on the depth 
Xi of the first interaction point. This is shown in Fig. Hlfor one hadronic interaction 
model. At high energy, and for air densities of relevance, almost all charged secondary 
pions always interact and all neutral pions decay immediately. Furthermore, the 
electromagnetic cascade initiated by the photons from vr" decay does not depend on 
the local air density. This makes the startup phase of a shower to a good approximation 
independent of Xi. 

The distribution Paxi depends, however, on the hadronic interaction model used 
for simulation. This is displayed in Fig. [3l both the mean AXi and the shape of 



(3) 
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Figure 4. Independence of PaXi of Xi, demonstrated for primary protons at 10^^ eV 
with QGSJetII. 



the distribution depend strongly on the chosen interaction model. In contrast, the 
distribution Pax2 exhibits only a small model dependence due to the universality of 
electromagnetic air showers after shower maximum, see Fig. [5] (left). 

An additional source of strong model dependence of Eq. ([3]) is the number of muons 
predicted for a given shower energy and depth. This can be seen in Fig. [5] (middle) where 
the frequency of finding a certain muon number for showers at 10^^ eV and different 
models is shown. The effect of this model dependence is displayed in Fig. [5] (right) by 
showing the folded distributions 

^effU^,;v, « / Pax, Pn, dE. (4) 

In real data analysis it must be taken into account that Xobs, -^e and iV^ are only 
known with a limited precision due to the detector and shower reconstruction resolution. 
In the model described by Eq. this would furthermore add the corresponding 
uncertainty distributions Pj^g and integrations over the three observables. The influence 
of the measurement resolution will be discussed below but is not included in ([3]) for sake 
of clarity. 



2.2. Optical telescopes 

Using fluorescence telescopes, the distribution of Xmax can be measured directly to 
obtain a handle on the value of the cross section at the highest energies [46, 47] . 

The direct observation of the position of the shower maximum allows us to simplify 
Eq. ([3]) by removing the term Pax, describing the shower development after the 
shower maximum and the distribution has to be replaced by a corresponding energy 
estimator based on the longitudinal shower proflle. The resulting distribution of Xmax 
can be written as 

Jl^^ /d«dEdX.dAX. ^^P^^,(E)P.^^ 

x*(X, + AX.-X„„)-^^j^^. (5) 
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Figure 5. Model dependence of longitudinal air shower development for proton 
primaries. Left panel: Slant depth distance from the shower maximum up to the 
depth where showers are attenuated to N^iE^ > 1 MeV) = 10^ electrons. Middle 
panel: Number of muons above 1 GeV after 1000 gcm^^ of shower development. Right 
panel: Same as left panel, but for showers selected according to 7.6 < logj^Q N^{E^ > 
1 GeV) < 7.7. While the plots on the left and in the middle are obtained from fixed 
energy simulations, the right plot is based on an energy spectrum oc E^^. 



Here E^m denotes the electromagnetic, i.e. calorimetric energy of the shower that 
can be obtained from the integration of the energy deposit profile. There is a very 
good correlation between the primary energy and the calorimetric energy [48]. This 
correlation and the weak energy dependence of the depth of shower maximum allow us 
to neglect Pe^^ = dpom/ d-Eem in our numerical studies. 

3. Analysis of cross section measurement methods 

The aim of all methods to measure the cross section is the determination of the 
interaction length Ajnt from measured distributions that represent the l.h.s. of Eqs. (!3|5|) . 

3.1. "k-factor" techniques 

The approximation of an exponential attenuation of the frequency of air showers after the 
penetration of large amounts of atmosphere is the basis of the fc-factor method [7,46,47], 
which has been used to analyze both Xmax and {Ne, Nfj) data. The exponential slope 
of the attenuation is typically denoted by A, which is then related to the hadronic 
interaction length by a so-called fc-factor 

A = k Ap_air • (6) 

While this method has been applied to data of air shower arrays as well as telescope 
detectors [7-12, 14], the definition of the fc-factor for a ground array, ks, and for a 
telescope detector, kx, are not identical. Air shower fluctuations enter differently into 
ks and kx and the detector resolutions are also very different. This can be understood 
based on Eq. ([3]) and ([5]), which can both be approximated by exponential distributions 
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Figure 6. Left panel: (iVo, 7V^)-method. Right panel: Xmax-tail method. The fitted 
open symbols show the dependence of the resulting exponential slope on the choice 
of the begin of the fit range. The arrow indicates the asymptotic behavior of the fit, 
deduced from Eq. (|lip . and its statistical uncertainty. 



for large Xmax, respectively Xobs- For the {Ne,Nf^) method this corresponds to 
diV^h 



dsec 9 

while for the Xmax-tail method it is 



OC e-^--/^- . (8) 

max 

This assumes that there are no significant non-exponential contributions to the tails 
of the distributions, which is not true in general as it will be shown in the following. 
Since the full distributions are described by Eq. ([3]) and ([5]) it is possible to calculate 
the tails based on the given convolution integrals. Each of the integrations will give a 
contribution, so for {Ne,Nfj^) one gets 

As = Aint ^AXi kAX2 k^cs = Aint with ks = kAXi kAX2 kres > (9) 

while for Xmax-tail just 

A-obs = Aint k^Xi ^res = Aint with kx = ^AXi kf^^ , (10) 

where Ajnt is the proton-air interaction length, k^Xi the contribution from the integration 

S /X 

of Paxi , kAX2 the one from Pax2 ^"^^ krL the part due to the detector resolution. The 
individual contributions of kAXi , kAX2 ^-nd krl^ are difficult to compute and not known 
in most of the existing analyses, except the one in Ref. [12] and even more complete in 
Ref. [14]. 

For a given experimental setup and analysis approach it is possible to estimate the 
values of kAXi, kAX2, ^res ^^nd kf^^ as well as their dependence on hadronic interaction 
models. This is demonstrated here for an artificial experimental setup as described 
below. For each of the hadronic interaction models QGSJetOIc, QGSJetII.3 and 
SIBYLL 2.1 a set of ~ 400,000 proton induced air showers is simulated with a 
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primary energy distribution (xE~^ between 18.7 < log^^Q-E/eV < 19.4 and a zenith 
angle distribution of dA^ oc cos ^^d cos 6*. Two histograms are produced and analyzed for 
every set of simulations: 

• Distribution of X^ax with 18.9 < log^Q-E/eV < 19.15. A Gaussian detector 
resolution with cr(Xmax) = 20gcm~^ and an energy resolution of a{E)/E = 
0.1 is assumed, which corresponds to values reported by the Pierre Auger 
Collaboration [49]. The resulting mean energy of this data sample is 10^^ ° eV. 

• Distribution of events with 1.26 < Ne/W^ < 3.98 and 0.28 < < 1.12 
at a vertical atmospheric depth of X^^l^ = 860gcm~^ (corresponding to the 
Pierre Auger Observatory) versus sec 6'. The A^e cut selects showers past their 
maximum (N^^^{10^^ eV) ~ 6.3-10^) as done, for example, in the pioneering Akeno 
analysis [7]% The detection resolution of A'e and A''^ are assumed to be Log- Normal 
with a resolution of (j(log]^o ^e)/ logio = 0.05 and cr(logiQ A"^)/ log^Q A"^ = 0.1 as 
well as the zenith angle uncertainty of a{6) = 1.0°. Due to significant model- 
differences in the prediction of muon numbers the selected data samples have 
slightly differing energies (from 10^^'^^ eV for QGSJetOI up to lO^^'^'^eV for 
SIBYLL). 

The histograms are fitted with an exponential function by a log-likelihood fit, which 
allows to correctly include empty bins in the analysis. The fit range is chosen as follows: 
While the end of the fit, X^^^, is always 50gcm~^ beyond the last non-zero entry in 
the histogram, the start of the fit, Xf^^^^, is varied. The resulting dependence of the 
exponential slopes on the start of the fitting range is then parameterized by 

A(Xf,i,J=Ao + Aie-^^?- , (11) 

where the asymptotic slope Aq for X^^^^^ ^ oo is taken to be the fit range- independent 
value of the exponential slope. In Fig. [6] this procedure is displayed for histograms 
obtained with the QGSJetII interaction model. It is found that the choice of the fit 
range has a major impact on the outcome of the slope. This is consistent with the results 
of [39]. For the (A'e,A'^)-method, generally no plateau is found and the resulting value 
for Aq is highly unstable. Thus, it is hardly possible to define a meaningful slope of the 
tail for this case. The situation for Xmax is somewhat better. While a stable plateau is 
still difficult to identify, Aq can be estimated reliably. With this analysis approach the 
resulting asymptotic values, Aq, are slightly lower than the values obtained for A directly 
from fits to the distributions. So neither the slopes nor the fc-factors can be directly 
compared to results from previous analyses. Nevertheless, the found model dependence 
and other methodical problems will apply in a similar way to previous analyses. The 
strong dependence of the slope on the chosen fitting range is a severe methodical problem 
of the fc-factor technique. Without a given method to infer a meaningful slope from 
such distributions, as for example the one proposed here, the fc-factor technique cannot 
produce reliable results. 

If In the recent analysis of the EAS-TOP Collaboration, showers are selected close to their maximum 
development [10]. 
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Table 1. Resulting fc-factors from CONEX study for primary protons at 10^^ eV. 
In each column the maximal and minimal fc-factors are highlighted. The difference 
between those maxima is denoted by A and is a measure of model-dependence. 



Model 


kAXi 


k\x2 


z,x 

"tcs 


lS 
"tcs 


fcx 


ks 


QgsjetII.3 
QgsjetOIc 

SiBYLL 2.1 


l.OOiO.Ol 
1.14±0.00 

1.02±0.00 


0.61±0.06 
0.48±0.05 
0.68±0.06 


0.92±0.01 
0.93±0.01 

0.93±0.01 


0.81±0.10 
0.50±0.14 

0.61±0.11 


0.92±0.01 
1.06±0.01 

0.94±0.01 


0.49±0.06 
0.27±0.08 

0.42±0.08 


Mean 
A 


1.05±0.00 
0.07±0.00 


0.59±0.03 
0.10±0.04 


0.92±0.00 
O.OliO.OO 


0.64±0.07 
0.16±0.09 


0.97±0.00 
0.07±0.00 


0.40±0.04 
0.11±0.05 



The obtained slopes Aq are then compared to the interaction mean free path length 
of the primary particle in the atmosphere, Ap_air, of the interaction model used in the 
CONEX simulations. In Tab. [1] the resulting /c-factors are listed together with their 
propagated statistical uncertainties. The total model induced uncertainties on the fcs/x- 
factors are ~ 7 % for the Xmax-tail and ~ 28 % for the {N^, Ai'^)-method. Looking to the 
fc-factor components one can see that a part of the model induced uncertainty enters 
already in the shower development up to the shower maximum (fcAXi), while the shower 
development after the shower maximum {k/^x2) adds the more significant amount of 
model-dependence. Another large contribution comes from /cf^g, which is mostly related 
to very differing model predictions on the number of muons. The factor /c^^ is only 
marginally model-dependent, which is originating from the very small model differences 
in the energy dependence of X^^^. 

While this study mostly serves illustrative purposes, and the chosen experimental 
and analysis parameters are arbitrary, it nevertheless demonstrates the impact of model 
dependence on fc-factor techniques. The found model dependence is rooted in the 
underlying air shower physics as well as in typical detector characteristics. 

3.2. Unfolding of the X^^y^- distribution 

An improvement of the cross section measurement techniques is achieved by explicitly 
accounting for air shower fluctuations [13]. This allows one to use not only the slope 
but also the shape of the Xmax-distribution. The measured Xmax-distribution, Eq. ([5]), 
is unfolded using a given Paxi -distribution to retrieve the original Xi-distribution. Of 
course, the AXi-distribution needs to be inferred from simulations, which ultimately 
introduces a comparable model dependence as in the fc-factor techniques [50]. The 
model dependence of the most important part of the kernel function, Paxi , can be seen 
in Fig. [31 

In the unfolding technique, a larger range of the Xmax-distribution is used. If the 
primary particles are known to be protons only, a cross section analysis can be done 
already with very limited shower statistics. On the other hand, the results of this method 
are more sensitive to a possible fraction of primary particles that are not protons. 
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4. An improved method to derive the proton-air cross section 

One of the shortcomings of cross section analysis methods applied so far is the missing 
connection between the cross section of the first interaction, that is to be measured, 
and the cross sections used in the calculation of the various probability distributions in 
Eqs. (ED and (El). 

In the following we consider the method of unfolding the X^ax distribution, 
applicable to fluorescence telescope measurements, and improve it by consistently 
accounting for the "measured" high-energy cross section. We will calculate the full 
shape of the distribution of -^max, depending only on properties of hadronic interactions 
at ultra-high energies. The impact of changing features of hadronic interactions at ultra- 
high energies on the longitudinal air shower development, i.e. P/\Xi , is parameterized and 
used within the calculation. It is straightforward to refine the modeling by incorporating 
detector acceptance as well as the energy distribution of analyzed data. 



4-1. Description of the X^^^^- distribution 

The description of the distribution of observed X^^x terms of o"p_air is based on 
Eq. ([5]), now written with the term for the experimental X^ax resolution 



oc / dXi dAXi dX^ax ^ PAXi(^,Ai„t,AXi + Xshift) 

Cl-^max J ^int 



X + AXi - X,„ax) PL{XZI. I ^max) • (12) 

The parameter Xshift allows us to shift the P^Xi -distribution by replacing AXi with 
AXi -|- Xghift- The introduction of the Xghift parameter is necessary in order to reduce 
the model dependence of the analysis [50]. This can be easily understood by looking 
at the PAXi-distributions shown in Fig. El which appear to be shifted between different 
models by up to ~ 60gcm^^. It is found that the Xshift parameter is highly model 
dependent and comprises many additional effects of the characteristics of high energy 
hadronic interactions on the Xmax-distribution. It can be demonstrated that differences 
between the inelasticity or the secondary multiplicity of the models may act as a cause for 
such an additional shift [51]. The only important assumption related to the role of Xghift 
in the cross section analysis is that any additional and unknown changing characteristics 
of hadronic interactions at extreme energies contributes mainly to a global shift of X^ax 
and thus AXi, while leaving the shape of the distributions unchanged. 

The cross section analysis proceeds then as follows. The X^^^-distribution is 
calculated from Eq. ( |T2l) for a given interaction length Aint and shift parameter Xshift 
and compared to the measured distribution. By performing a log-likelihood fit with the 
interaction length and an overall shift in depth as free parameters, the cross section and 
the la uncertainty band is found [50]. 

One major improvement with respect to previous cross section analysis approaches 
is the consideration of the impact of a changing cross section on the resulting shower 
development described by Paxi • It is assumed that the cross section is reasonably well 
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Figure 7. Modified extrapolation of the proton-air cross section. Shown are the 
original prediction from SIBYLL together with a 20 % increased respectively decreased 
cross section at 10"'^^ eV. 



known at the energy corresponding to that of the Tevatron colhder. Starting from this 
energy of ~ 2 x 10^^ eV, cross sections of the model used for the calculation of Paxi are 
rescaled by a factor f{E) that increases logarithmically with energy and ensures that 
the modified model cross section matches the corresponding value of the fit parameter 
Aint at the considered shower energy. Here we will use 10^^ eV as reference for the scaling 
parameter being fig at this energy. Then the scaling factor reads 

fiE)-l + if,, - 1) 1^(1019 eV/lOi^eV)' ^^^^ 

for E > 10^^ eV and f{E) = 1 otherwise. The idea of the rescaling of the cross section is 
shown in Fig. U\ where the proton-air cross section of SIBYLL is shown together with 
cross sections that are scaled up and down by 20 % at 10^^ eV. 

So far the dependence of Paxi on the cross section has been neglected in air shower 
based cross section measurement. In fact, any attempt to determine the proton- air cross 
section without taking this dependence into account is overestimating the impact of 
o"p-air on the analyzed observables since part of the measured effect must be attributed 
to a modified development of the air shower and not to the fluctuations of the first 
interaction point. 

To implement the idea of modified cross sections in the data analysis one needs 
to parameterize the AXi -distribution and determine its energy dependence and the 
modification of this distribution as function of the cross section scaling parameter. This 
is done using a Moyal distribution extended by one parameter and described in detail 
in Appendix A In Fig. [8] example AXi-distributions are shown together with fits of the 



extended Moyal distribution and also the final parameterizations. 
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AX^ [gcm=] AX, [gem"'] 

(a) QGSJetII.3 (b) SIBYLL 2.1 



Figure 8. Simulated AXi-distributions at 10^^ eV for two hadronic interaction 
models. Included are the fits of the extended Moyal distribution Eq. (jA.ip as well 
as the final parameterization. 

4-2. Fitting range and stability 

To investigate the sensitivity of the method to the choice of the X^ax-range used for the 
fit, we generate sets of 3000 simulated data showers and analyzes them with Eq. ( I12p . 
The Gaussian detector resolution for X^ax is chosen to be 20gcm~^, which is the value 
reported by the Pierre Auger Collaboration [49]. 

The interval of the Xmax-distribution used to fit the model is important for the 
resulting statistical as well as systematical uncertainties. The first aspect is obvious 
since a reduction of the dataset is clearly leading to a reduced statistical power of the 
reconstruction. The latter one is mostly because of the possible contamination of the 
-^max-distribution by cosmic ray primaries other than protons. All primary nuclei heavier 
than protons are producing shallower X^ax compared to protons. Primary photons, on 
the other hand, are deeply penetrating and have a larger X^ax than protons. A restricted 
fitting range in X^ax can thus be used to enrich the considered fraction of protons and 
reduce a possible contamination by other cosmic ray primaries. 

The resulting impact of the chosen fitting range on the performance of the 
reconstruction is shown in Figure [H The position of the peak of the Xmax-distribution, 
Xpeak, is used as a reference to define the fitting range. The starting point as well as 
the ending point of the fit are expressed only relative to Xpeak- 

Evidently, it is beneficial to chose the fitting range as large as possible to get the 
smallest resulting statistical uncertainties. It is possible to shift the beginning of the 
fitting range relatively close to Xpeak without loosing too much statistical power (cf. 
Figure [9] (a)). It is found that, with the used parameterizations, the beginning of the 
fitting range can be set to 50 gcm~^ in front of Xpeak- In the following this is the adopted 
default choice. 

The choice of the end of the fitting range has a similar impact on the reconstruction. 
In Figure [9] (b) it is shown how the reconstruction is degrading, while choosing a shorter 
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(a) Start of the fitting range, ^max* (relative to the peak of the Xmax-distribution) 
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(b) End of the fitting range, X^°^ (relative to the peak of the X,nax-distribution) 



Figure 9. Impact of the chosen fitting range in X^^^ on the resulting cross section 
(left) and Xshut (right) at 10^^ eV. Each point corresponds to the mean of 100 
reconstructions and the error bars denote the resulting RMS. Both, the reconstructed 
-'^max as well as the AXi-distribution are produced independently with the specified 
interaction model. The lines are just to guide the eye. 



fitting range. Since the photon fraction of ultra-high energy cosmic ray primaries is 
already strongly constrained [52] there is no special need to restrict the upper end of 
the fitting range. In the following the upper end of the fitting range for the log-likelihood 
fit of the Xmax-distribution is set to the value of the maximum Xmax plus 40gcm~^. 

The cross section dependent parameterizations of the AXi-distributions are 
introducing a slight systematic overestimation of the reconstructed cross sections of 
10 — 20 mb, corresponding to < 5%. This can be taken into account within the 
determination of the total systematic uncertainties of the measurement. The statistical 
resolution of the reconstruction is around 20 mb for 3000 events. 

5. Comparison of the performance of different analysis methods 



In order to test the ability of a reconstruction method to recover a changing input cross 
section, air shower simulations with a modified cross section (cf. Appendix A. 2 ) are 
performed. The simulated air shower data are then reconstructed and the found cross 
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Figure 10. Resulting systematics and dependence on models of cross section 
reconstruction techniques. Every analysis was performed 100 times on 3000 simulated 
EAS; Shown are the mean results. The central shaded area shows the ability of a 
method to reconstruct the input cross section of simulated air shower data with the 
same model that was used for the simulations. The thick line marks the mean result for 
all used models, which are QGSJetOIc, QGSJetII.3 and SIBYLL2.1. The width of 
the central band is caused by statistical fluctuations. The outer shaded area denotes 
the maximum deviation of the cross section reconstructed with a model that differs 
from the one used to generated the simulated dataset. 
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Figure 11. Shower-to-shower fluctuations of X^ax (left) and log^Q A'^o (right) for 
primary protons. The solid line is the original model prediction of QGSJetII, while 
the dashed line is for showers having their first interaction at the fixed slant depth 
corresponding to the predicted proton-air interaction length of QGSJetII. The dotted 
line is for simulated air showers with increased cross section by a factor of two, /ig = 2, 
in addition to the fixed first interaction point. 



sections (Tree are compared to the modified input cross sections (Tmodified = /i9 Cmodei- 
The result of this analysis for all cross section reconstruction methods is summarized in 
Fig. [TOl The central shaded area together with the solid line demonstrates the principal 
self-consistency of the analysis technique. The outer shaded area can be interpreted as 
the maximum possible mo del- dependence of the reconstruction method. 

With the (A'^e, A''^)-method it is very difficult to reconstruct a modified cross section. 
Additionally the model dependence of the result is generally large. It is difficult to 
quantify, but can easily be of the order of hundreds of mb. 

A very much better performance is achieved by the Xmax-tail analysis. The 
deviation of the reconstruction cross section from the input value is not getting larger 
than 50 mb and the model-dependence is between 100 — 200 mb. 

The Xinax-unfolding technique can find cross sections that are not deviating 
much from the original model value. However, it shows a clear systematic trend of 
underestimating small cross sections and strongly overestimating large cross sections. 
Also the model dependence is not negligible, ranging from ~ 75 mb at very low cross 
sections, ~ 200 mb at intermediate cross sections to many hundreds of mb at large cross 
sections. 

Finally, the Xmax-model demonstrates its unique ability to retrieve very consistently 
modified cross sections over a wide range. The model dependence is smaller than for 
the Xjnax-tail method and is < 50 mb at low cross section while it grows to ~ 100 mb at 
large cross sections. 

The generally better sensitivity to reconstruct smaller compared to larger cross 
section is resulting from the increasing importance of the fiuctuations of Xi for the final 
Xmax-distribution. This facilitates the measurement of a small cross section, while for 
large cross sections Xmax-distributions are mostly shaped by fiuctuations of the shower 
development Paxi and the detector resolution, making a measurement of the cross 
section generally more difficult. 
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Due to the combined impact of the fluctuation of both electron as well as muon 
numbers, the remaining sensitivity of the [Ng, iV^)-method to the proton-air cross section 
is very much reduced compared to methods using the observable Xmax, see Fig. [HI 

It is not straightforward to quantify the systematic effect caused by interaction 
models on the results of the proton-air cross section analyses. While the magnitude of 
the outer uncertainty band in Fig. [10] suggests relatively large model dependence of the 
results, this just reflects the existing incompatibility of the underlying models. Thus, it 
cannot be used to make any conclusion on the comparison of real data to the different 
models. For the analysis of real data, model dependence should be evaluated by the 
independent reconstruction of the data based on all available interaction models, and 
subsequent quoting of the mean result. The interval spanned by the reconstructions 
can then be used as an estimation of the systematic uncertainty induced by the models. 
This results in a smaller dependence on interaction models compared to Fig. [THl 

6. Summary 

It is demonstrated in methodical studies that all techniques to determine the proton-air 
cross section from air shower data are in fact based on the same fundamental formulation 
that describes the longitudinal air shower development. The requirement of using air 
shower Monte Carlo simulations to interpret EAS observables inevitably introduces a 
dependence on the hadronic interaction models used for the simulations. 

It is found that the magnitude of the model dependence is similar for all cross 
section reconstruction methods; Only the {Ne, iV^) -technique exhibits a significantly 
larger model dependence, which is due to the additional strong model dependence on 
the prediction of muon numbers. In addition, since the prediction of electron and muon 
numbers are both depending on interaction models, it is not possible to define model 
independent cuts for the selection of air shower events. Thus, already on the level of 
event selection, a model dependence is introduced in {Ne, iV^)-analyses. 

All fc-factor based techniques are suffering from non-exponential contributions to 
the slope of the Xmax- respectively {Ne, Nfj) frequency-distribution at large depths. Thus 
the exponential approximation used to define fc-factors, is only of limited accuracy. Any 
non-exponential contribution creates a strong dependence of the exponential slope on 
the chosen fitting range (see also Ref. [53]). 

Generally, fc-factors are depending on the resolution of the experiment and can 
therefore not be simply transferred from one experiment to another, as done in some 
recent analysis [54]. In particular, fcx-factors are inherently different from /cs-factors 
and can therefore not be transferred from an Xmax-tail analysis to that of ground based 
frequency attenuation or vice versa. 

An Xjnax-based analysis technique, as proposed in Ref. [13], was further improved 
to reduce methodical and model induced systematic uncertainties on the cross section 
result as much as possible. The presented improved method is a model of the X^ax- 
distribution based on simple considerations on longitudinal air shower development 
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combined with the parameterized impact of a changing cross section on the resulting 
air shower development. Previous reconstruction techniques assume a static AXi- 
distribution or fc-factor, independent of the modified cross section. As, in general, the 
reconstructed cross section deviates from the predicted value of the hadronic interaction 
model used to generate the AXi-distribution or fc-factor, a discrepancy of the data with 
the air shower simulations is inevitable. The data can not be described by the models 
used for the air shower simulations and it is not sufficient to adopt the measured cross 
section just for the first interaction point, while keeping the rest of the shower simulations 
unchanged. 

For the case of a proton-dominated cosmic ray data sample around 10^^ eV, which 
was studied in this article, it is possible to determine the proton-air cross section 
with good precision. The statistical uncertainties of the analysis of samples with 
3000 events are about 20 mb (5 %) and the systematic uncertainties introduced by the 
parameterizations of the Xmax-inodel are of the same order of magnitude. The model- 
induced systematic uncertainty is between 50 and 100 mb, corresponding to relative 
uncertainties of < 10 % to ~ 20 %. A further reduction of this model-related uncertainty 
is possible, but needs improvements in the understanding of the underlying hadronic 
interaction physics. 

Furthermore, the second free parameter of the new analysis technique, Xghift, is 
a further handle that is sensitive to the physics of hadronic interactions at ultra-high 
energies. A non-zero result of Xsmt indicates the deficiency of the underlying hadronic 
interaction model to describe the measured data distribution. It can be interpreted 
in terms of a modified characteristics of secondary particle production in hadronic 
interactions [51]. 
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A modified version of the Moyal distribution is introduced to parameterize the AXi- 
distribution 



Pi(AXi) = -— =e V J with z = )- . (A.l) 

p y/ Ztt P 

The additional third degree of freedom 7 makes the parameterization more flexible, 

allowing us to improve the description around the peak of the AXi-distributions. The 

normalization value is not a free parameter, but is numerically computed to correctly 

normalize the distribution. Since the Moyal function with the two parameters a and f3 

itself is normalized, we only have to correct for the 7 parameter 

1 

0.965 + 3.685 x 10-3 (7 + 0.366)" 



Appendix A.l. Energy dependence 

To determine the energy dependence of the parameters a, j3 and 7, the parameterization 
Eq. (lA.ll) is fitted to AXi-distributions generated by CONEX at several primary 
energies. The energy dependence can then be interpolated by a polynomial of 2°*^ 
degree 

a\E) =al + al \og^^{E/eY) + al \ogl^{E/eV) 

f5\E) = PI + PI \og,,{E/eV) + PI \ogl,{E/eV) (A.3) 
f{E) = 7o^ + 7^ \og,,{E/eY) + 72^ \ogl,{E/eY) . 

The chosen polynomial interpolation is well suited to reproduce the found 
dependences over the considered energy range. The results are illustrated in Fig. lAll 
and the parameters are listed in Table lAll A first result is that even the energy 
dependence of AAi is depending on the high energy interaction model used during 
air shower simulations. 



Table Al. Parameters for describing the energy dependence of AXi-distributions. 



Model 


Index i 


al [g/cm] 


/3| [g/cm] 


7| [1] 







-871.25±27.76 


85.15±14.49 


5.94±1.43 


QGSJetOIc 


1 


113.91±3.01 


-7.49±1.57 


-0.71±0.16 




2 


-1.64±0.08 


0.22±0.04 


0.02±0.00 







-1320.67±24.86 


226.24±11.68 


6.04±0.39 


QGSJetII.3 


1 


165.29±2.70 


-22.31±1.27 


-0.63±0.04 




2 


-3.04±0.07 


0.59±0.03 


0.02±0.00 


SIBYLL 2.1 





-606.32±20.37 


131.27±9.56 


0.33±0.43 




1 


79.17±2.20 


-11.78±1.03 


-0.05±0.05 




2 


-0.45±0.06 


0.30±0.03 


O.OOiO.OO 
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Figure Al. Parameterization of tlie energy dependence of the AXi-distributions. 
The lines denote the polynomial interpolations given in the text. 



Appendix A. 2. Cross section dependence 

To infer the dependence of Paxi on Ap_air, a modified CONEX version [51] was used 
to generate AXi-distributions for several values of fig, see Eq. ( 1T3|) . The energy 
dependence of all hadronic cross sections, i.e. p-, vr-, and K-air interactions are then 
altered by the modification factor. This implies that the model uncertainties above 
10^^ eV are rising proportional to the logarithm of the energy. This choice seems 
certainly reasonable, however, in the future also other energy dependences of f{E) 
may turn out to be useful. 

The resulting dependence of AXi on a changing cross section for air showers at 
10^^ eV is shown in Fig. IA2[ 

The dependence of the a, (3 and 7 parameters of the modified Moyal distribution 
on /i9 is parameterized as follows 

a'^(/i9)=a^ + ar/i9 + a2/i9 

/5'^(/l9)=/?o+/?l7(/l9-/?2) (A.4) 



Table A2. Parameters for describing the cross section dependence of AXi- 
distributions. 



Model 


Index i 


a.1 [g/cm] 


/3? [g/cm] 


li [1] 







730.09±0.25 


10.33±0.00 


0.03±0.03 


QGSJetOIc 


1 


-33.61±0.34 


11.83±0.00 


0.10±0.03 




2 


5.07±0.09 


-0.13±0.00 


0.20±0.07 




3 






-O.OliO.Ol 


QGSJetII.3 





754.10±0.27 


8.46±0.08 


-0.25±0.02 




1 


-35.86±0.35 


8.28±0.19 


0.10±0.02 




2 


5.14±0.09 


-0.08±0.01 


-0.20±0.09 




3 






0.02±0.00 


SIBYLL 2.1 





774.36±0.27 


11.67±0.09 


-0.18±0.04 




1 


-46.36±0.36 


4.76±0.16 


0.10±0.04 




2 


7.12±0.10 


0.12±0.01 


-0.00±0.15 




3 






0.04±0.01 
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■lOEeV 'lOEeV 'lOEeV 

Figure A2. Parameterization of tlie /19-dependence of tlie AXi-parameterizations 
witli respect to a changing cross section. The lines denote the interpolations as 
discussed in the text. 



r{fl9) = 7o + 7l7(/l9 - 72 ) + 73 /l9 ■ 

The resulting parameters of the interpolations are listed in Table IA2[ Furthermore, it 
is convenient to introduce the relative change as 

A«'^(/i9)=a'^(/i9)-«'^(l) 

A/5'^(/i9)=/3'^(/i9)-r(l) (A.5) 
A7'^(/i9) =7'^(/i9)-7'^(l) • 

Appendix A. 3. Parameterization of the AX i- distribution 

Combining the found dependence of AXi-distributions on the energy E and the cross 
section modifier /19, it is now possible to construct a parameterization of Paxi in terms of 
E and cr{E). Firstly, the principle energy dependence of Paxi is evaluated by calculating 
a'^{E), P'^{E) and 7^(-E) using Eq. flA.3l) and the results from Tab. lAll Secondly, the 
effect of the changed cross section is added. To achieve this, the first step is to calculate 
the corresponding deviation of the cross section a^E) extrapolated to 10^^ eV assuming 
Eq. (fT3l) . The resulting factor fig is then used to evaluate Aa°"(/i9), ^/^'^(/ig) and 
A7'^(/i9) using Eq. flA.51) and the results listed in Tab. IA2[ In accordance to Eq. f|T3l) 
the energy dependence is assumed to be logarithmic and vanishing below 10^^ eV. This 
yields the final set of parameters of Paxi 

a{E, a{E)) = a%E) + F{E) Aa'^if^g) 

(3{E, a{E)) = (3\E) + F{E) A[3^{fig) (A.6) 
j{E,aiE))=^%E) + A77/19), 

with 

F(E) = l^ ^<10^'eV 
^ ' ^ ln(E/lPeV)/ln(10EeV/lPeV) E > lO^^eV ' ^ 
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